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ABSTRACT 

We investigate properties of plasma turbulence from magneto-hydro dynamic (MHD) to sub-ion 
scales by means of two-dimensional, high-resolution hybrid particle-in-cell simulations. We impose an 
initial ambient magnetic field, perpendicular to the simulation box, and we add a spectrum of large- 
scale magnetic and kinetic fluctuations, with energy equipartition and vanishing correlation. Once the 
turbulence is fully developed, we observe a MHD inertial range, where the spectra of the perpendicular 
magnetic field and the perpendicular proton bulk velocity fluctuations exhibit power-law scaling with 
spectral indices of —5/3 and —3/2, respectively. This behavior is extended over a full decade in 
wavevectors and is very stable in time. A transition is observed around proton scales. At sub-ion 
scales, both spectra steepen, with the former still following a power law with a spectral index of ^ —3. 

A —2.8 slope is observed in the density and parallel magnetic fluctuations, highlighting the presence of 
compressive effects at kinetic scales. The spectrum of the perpendicular electric fluctuations follows 
that of the proton bulk velocity at MHD scales, and flattens at small scales. All these features, 
which we carefully tested against variations of many parameters, are in good agreement with solar 
wind observations. The turbulent cascade leads to on overall proton energization with similar heating 
rates in the parallel and perpendicular directions. While the parallel proton heating is found to 
be independent on the resistivity, the number of particles per cell and the resolution employed, the 
perpendicular proton temperature strongly depends on these parameters. 

Subject headings: plasmas - solar wind - turbulence 


1. INTRODUCTION 

Turbulence is an ubiquitous phenomenon in space and 
astrophysical plasmas. Although it is generally driven 
by violent events or instabilities at large scales, a further 
cascade is responsible for transferring energy via non¬ 
linear coupling from the large injection scale to much 
smaller scales, through the ion and the electron charac¬ 
teristic regimes, where they are eventually dissipated. In- 
situ measurements in the solar wind represent a unique 
opportunity to study those processes, since they pro¬ 
vide observations in a huge range of sca l es (see for ex¬ 


ampl e reviews by Tn fc Marsch 
Velli ( 2011 ); Alexandrova et al. 
bonel (2013)). The estimated tur 


.g., iMacBride et ^ 2008 


rate _ 

Hellinger et al. ' 201l[ 2013 


1995 ); Matthaeus fc 
013p ; Bruno Car 


Dulent energy cascade 


Cranmer et al. 2009 
is comparable to the pro- 
ton heating needed to explain the non adiabatic evolu- 
tion of the solar wind plasma during its expansion (e.g., 
Marsch et al.||2004| |Matteini et al.|[2007|). This suggests 
that turbulence plays an active role m transferring energy 
from electromagnetic fields to particles and heats the so¬ 
lar wind plasma. However, the processes that ultimately 
lead to heating in a collisionless turbulent medium are 
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still unknown. 

Supporting evidence of a turbulent cascade is provided 
by the observed energy spectra, which exhibit a power- 
law behavior over a large range of scales, spanning nearly 
four decades in frequency. The spectral index of mag¬ 
netic and kinetic spectra varies with the temperatur e 


of the solar wind streams (Grappin et al. 1990 1991) 
The latter is in turn correlated with the stream speed 
and, although to a smaller extent, with the degree of 
Alfvenicity, i.e., th e correlation between kinetic and mag¬ 
netic fl uctuations ( Podesta fc Borovsky|2QlQ Chen et ^ 
2013a). However, on average, at tluid-like scales a typ¬ 
ical Kolmogorov power law with spectral index —5/3 is 
usually observed for magnetic fluctuations, while kinetic 


ing (Podesta et al.|2006l 2007 

iTessein et al.|2009 

Salem 

et al.|:^0u0| 

; kJhen etal.EOTTFr 

Wicks et alFAllk. ] 

.n par- 


balanced turbu lence, i.e., zero cross-helicity (Podesta & 


Borovsky 2010). In the same range of scales, a certain 
amount of residual energy, i.e., an excess of magnetic to 
kinetic energy, is typically observed, followin g a well de¬ 
fined p ower-law scaling with an index of —2 (Chen et al. 
2013a|). The electric field spectrum is observed to follow 
the vel ocity spectrum, wh en measured in the solar wind 
frame (Chen et al.||2011a), while density fluctuations ex¬ 
hibit a Kolmogorov-like cascade. 

In the vicinity of the ion inertial length scale, a break 
in the magnetic field power spectrum is observ e d (|Bein- 
roth fc Nenbauer 1981| |Goldstein et al.| 1994[ Leamon 
et al. 1998| 1999p. Early observations of the spectrum 
of magnetic huctnations in a restricted region above the 
break found a power-law scaling wit h a variable spectral 
index, ranging from —2 to —4 (e.g., Leamon et al.||199'8' 
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1999| 

Bale et al.||2005t |Smith et al.||2006t Alexandrova] 

et al.| 

2008b|al Kiyani et al.j 

2009] jSahraoui et al. 2009; 

Chen et akpOlOj Salem et a 

.[[2012p. However, more re- 


cently, observations extended to smaller scales suggest 
a general conve rgence of the spectr a towards a spectral 
index of —2.8 dKiyani et al. 2009 Alexandrova et al 


2009 Sahraoui et al. pT^ ) , or towards a power-law scal¬ 


i ng ot —8/3, exponentially damped at sub-electron scales 


( [Alexandrova et a^ | 2012[ ). Magnetic fluctuations at sub 
ali 


proton scales are also cnaracteriz ed by a reduction of 


the m agnetic variance anisotropy (Podesta & TenBarge 
2012|), and by an increase of the niagnetic compressibil¬ 


ity ([Alexandrova et al.||2008a| [Salem et al.||2012t [Kiyani 

gtll.[[2013| ), suggesting a change in the nonlinear mter- 
actions ruling the cascade. This is partially confirmed by 
the measured increase of the intermittency at ion scales 
([Alexandrova et al.[[20Q8a1 [K iyani et al.[[2009[ [2013| [Wu 
et al.[|2013 ' (Jhen et al.|[2014 ), although a clear behavior 
ot the flatness at smaller, sub-ion, scales has not been 
identified yet. 

There are observational indications that, below the ion 
inertial length scale, the electric field spectrum decou¬ 
ples from the velo city field and flattens ( Bale et al.|2QQ5 


Salem et al. 2012) but, due to the high noise level, present 
data do not allow to determine the existence of a power- 
law scaling at sub-ion scales. Density fluctuations show 
a plateau just before the ion scales, while they follow a 
power law between the ion and the electron scales, with 
the same s pectral index as the one of the magnetic field 


spectrum (Chen et al. 


Properties of turbu. 


2012 2013b) 


ence have been extensively ana¬ 


lyzed by means of direct numerical simulations (DNS), 
employing many different methods and models. Al¬ 
though several features of the solar wind turbulence can 
be partially recovered, we are still far from a compre¬ 
hensive picture. At large fluid-like scales, DNS of in¬ 
compressible MHD and reduced MHD (RMHD) return 


- ^ -- --- ---- 

or —3/2 (e.g., Maron & Goldreich 2001[ 

i [Muller et al. 

2003] jlVIuller Grappin||2UU5[ IVIason et < 

al. |2UU8| jPerez 

& Boldyrev||2009| jBeresnyak &; Lazarian 2009] Grappin 

& Muller|20101|Boldyrev et al.|20111|Lee et 

al. 12010] Chen 

et al. 201 lb| Beresnyak 2011] Perez et al. 

|2012p. These 


of the nonlinear interactions regulating the cascade and 
the cascade rate. Moreover, within the inertial range, a 


transition between different regimes can occ ur (Mininni 


& Pouquet [2007 Verdini & Grappin 2012). More so¬ 
phisticated DN S, i ncluding other phy sical processes like 
expansion effects ([Dong et al.[[2C)14[), Hall-MHD (e.g 

^ i ;u _ 1 ' ii ll onno L ’ 


Matthaeus et al.| 

Shukla||2009| 




([Gomez et al.|20l3), gyrokinetic 


nr 


20031 IGomez et ai] 2008| [Shaikh fc 
ank||20 09[), reduced Hall MHD 
Hx 


lowes 


et al.|200^ Ten- 


Barge et al. 1201311, and hybrid particle -in-cell (PIC) sim- 


ulations (jVasquez & IVlarkovskii 2012), all produce spec 
tral indices consistent with —5/3. Anyway, the restricted 
width of the inertial range prevents firm conclusions. 

As far as the small kinetic scales are concerned, DNS 
including proton and electron physics return a qualita¬ 
tively unified picture. At sub-proton scales, they re¬ 
produce an increase of the ratio of the electric to mag¬ 
netic power, togeth er with a flattening of the electric 
field spectr um (e .g., Dmitruk fc Matthaeus[[2006 ; Howes 


et al. 120081 |20Tli |; ^ervidio et al. | :^0l:^l [ G6mezet al. | :^0l;^ 


Perrone et ah] 2013] ]Parashar et al.]]20141 

Passot et ah] 

2014] Valentini et al. 2014] Servidio et al. 

2015p, and a 


power near the ion scales (e.g., Matthaeus et al. 

2003 

Dmitruk & Matthaeus 2006] ]Parashar et al. 


]Ser- 

vidio et al. 12012] ]Vasquez & Markovskii|20121 

]Rodriguez 

Imazio et al. 2013] Valentini et al. 2014p. ] 

However, a 


Hall MHD ( 

Shaikh & Shukla 

2009 

Martin et al. 2013]), 

Electron-MHD (]Biskamp et a. 


Ngetal.|2U03||(Jho| 

& Lazarian 

|2004l |2UUy Shaikh 2Ue 

)9), and gyrokinetic 

(]Howes et a 

I. 2008^ reported 

a spectral index of —7/3 


steeper spectra have als o been observed: a spectral index 
of —2.8 in gyrokinetic ([Howes et gil. 2Qllj 


TenBarge & 


Howes 1 2 (H 3 TenBarge et al.|2U13p a nd hmteXarrnor 


dius (FLR)-Landau fluid simulations (Passot et al.|2Q14 ) 
or a —8/3 power law both in 3D electron-IVIHD ([iVleyrand 

—- fa - •- --- -- — 


& Galtierl 


an d in s trong kinetic-Alfven turbulence 


([Boldyrev Perez[[2U12[). Magnetic spectral indices in- 
between about —2.6 and —3 have also been observed in 


full PIG simulation s (e.g., [Gamporeale fc Burgess[[2Qll 


Ghang et al.||2Qlll |Wan et al.||2U12[ [Karimabadi et al. 


20131 |Wu et al.[[20i3p. 

In situ measurements of the proton velocity distribu¬ 
tion functions show the presence of an ubiquitous tem¬ 
perature anisotropy between the directio n parallel and 


perpendicular to the mean magnetic field (Marsch et al 
1982 Hellinger et al. 2006), and a non adiabatic e volution 


thus suggesting, as 
already mentioned, an active role played by the turbu- 


of theTolar wind plasma during i ts expansion (Marsch 
et al.l 2004| [Matteini et al.[[2007[) 

id' 


lence in exchanging energy between fields and particles. 
Hybrid PIC simulations have shown an overall (macro¬ 
scopic) collisionless proton heating, with signatures of a 
preferential proton heating in the perpendicular direction 


Parashar et al. 

]2009| ]Markovskii et al.]|2010| 

Markovskii 

& Vasquez 

201] 

Ll]Vasquez & Markovskii|20121 

Verscharen 

et al. 2012 

Parashar et al. 2014|). Vlasov-H 

ybrid simu- 


as temperature anisotropies, can be produced by the tur¬ 
bulence, mostly concentrated in re gions near and around 


the peaks of the current density ([Servidio et al. 2012 
Valentini et al.[[2014| [Perrone et al.||2014| |bervidio et al'. 

mT^'. 

Inour previous w ork (Franci et al.|2015) (named here¬ 
after as Letter 1[) , we presented results from a high- 
resolution hybrid (fluid electrons, PIC protons) two- 
dimensional (2D) simulations of turbulence. The spectra 
of various fluctuations (magnetic, kinetic, density and 
electric field), along with the magnetic compressibility 
and the non-dimensional ratio of the density and the 
magnetic fluctuations, simultaneously matched several 
features observed in the solar wind. In particular, for 
the magnetic field we showed that high-resolution hy¬ 
brid simulations, although limited to a 2D geometry, are 
able to capture the nonlinear dynamics at fluid-like MHD 
scales and at subproton scales, both within the same nu¬ 
merical domain. In this paper, we analyze in further 
detail the spectral properties of several fields, also show¬ 
ing their stability with time. Moreover, we investigate 
the shape of the electric field spectrum, by estimating 
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the separated contributions from different terms in the 
generalized Ohm’s law. Finally, we study the proton tem¬ 
perature anisotropy and the proton heating, also quanti¬ 
fying the dependence from the resistivity coefficient and 
the number of particles-per-cell (ppc) employed in the 
simulations. 

The paper is organized as follows: In Section we 
describe the numerical setup employed, define the phys¬ 
ical units and normalizations in the code, and provide 
the parameters of our initial conditions. In Section 
we describe the results of the performed simulations. In 
Section]^ we validate such results, by investigating the 
importance of a careful choice of some relevant numeri¬ 
cal parameters. Finally, in Section we summarize the 
achievements of our simulations and^discuss them in the 
framework of both observational and previous numerical 
and theoretical studies. 


2 . NUMERICAL SETUP AND INITIAL 
CONDITIONS 


We make use of a 2 D hybrid code, where electrons are 
considered as a massless, charge neutralizing fluid with 
a constant temperature, whereas ions are described by a 
PIC model and are advanced by the Boris’ scheme, which 
requires the fields to be known at a half time step ahead 
of the particle velocities. This is achieved by advancing 
the current density to this time step with only one com- 
pntat ional pass throu gh the particle data at each time 
step ( Mattlie^|1994 ). 

The characteristic spatial and temporal units used in 
this model are the proton inertial length dp = cjup^ Up = 
( 47 rne^/mp)^/^ being the proton plasma frequency, and 
the inverse proton gyrofrequency = {eB^^jmpC)~^^ 
respectively. Magnetic fields are expressed in units of the 
magnitude of the ambient magnetic field, i.e., 5o, while 
velocities are expressed in units of the Alfven velocity. 


i.e., va = 5o/(47rnmp)^/^. The plasma beta for a given 
plasma species, protons (p) or electrons (e), is Pp^e = 
SirnKBTp^e/B q. Quantities and symbols used in these 
definitions are: the speed of light, c, the number density, 
n, which is assumed to be equal for proton and electrons 
{up = He = n), the magnitude of the electronic charge, 
e, the proton mass, m^, the Boltzmann’s constant, Kb, 
and the proton and electron temperatures, Tp^^. 

The 2D computational domain lies in the (x, y) plane, 
while the ambient magnetic field Bq is along the 2 :- 
direction. Accordingly, each field T will be decomposed 
in its perpendicular (in-plane) component, and its 
parallel (out-of-plane, along z) component, Ty, with re¬ 
spect to Bq. The only exceptions will be the proton beta 
and temperature, for which T and || will denote direc¬ 
tions with respect to the local magnetic field. 

The adopted simulation box is a 2048^ square grid. We 
tested different resolutions {Ax = Ay = 0.5, 0.25 and 
0.125dp), and consequently different box sizes (Lbox = 
1024, 512 and 256dp), as well as different numbers of 
ppc, ranging from 500 to 8000 (see Table [^. The time 
step for the particle advance is At = 0.025 whereas 
the magnetic field B is advanced with a smaller time 
step. Ate = At/10. 

Initially, we assume a uniform number density n = 1, 
a proton parallel beta /3p\\ = 0.5, and a temperature 
anisotropy Ap = Tp±/Tp\\ = 1. Electrons are isotropic. 


Run 

Ax 

(dp) 

Abox 

(dp) 

V 

(47r/wp) 

ppc 

A 

0.125 

256 

5 X 10-* 

8000 

B 

0.125 

256 

5 X 10“* 

4000 

C 

0.125 

256 

5 X 10-* 

2000 

D 

0.125 

256 

5 X 10-4 

1000 

E 

0.125 

256 

5 X 10-4 

500 

F 

0.125 

256 

1 X 10-4 

8000 

G 

0.125 

256 

1 X 10-3 

8000 

H 

0.250 

512 

1 X 10-3 

8000 

I 

0.500 

1024 

2 X 10-3 

8000 


TABLE 1 

List of simulations and their relevant parameters 


with Pe = 0.5. We add an initial spectrum of magnetic 
and velocity fluctuations in the {x^y) plane, composed 
of modes with — 0.2 < kx^y < 0.2 in each direction and 
random phases. These initial fluctuations are charac¬ 
terized by energy equipartition and vanishing correla¬ 
tion between kinetic and magnetic fluctuations, and their 
global amplitude is B^^^ ^ 0.24. Hereafter, 

^rms ^ 

will denote the root mean square value (rms) of a quan¬ 
tity T, with (T) being its space-averaged value over the 
whole 2D simulation domain. The initial magnetic fluc¬ 
tuations can be expressed in the form 

-^t(t, y ) — ~ ^ ^ [-^t(^cc: ky) 

^ K,ky ( 2 ) 

X ex.p{i{kxX + kyy + ^{kx, ky))) -h c.c.]. 


The initial bulk velocity fluctuations u_i{t = 0) are as¬ 
sumed to have the same form as in Eq. Q , with different 
random phases. 

We introduce two dimensionless quantities, i.e., the 
normalized cross helicity, ac, and the normalized resid¬ 
ual energy, an: 


crc{x,y) 

crR{x,y) 


2 u B 

WTW 


( 3 ) 

( 4 ) 


which define the two-dimensional geometry of the fluc¬ 
tuations. With the initial conditions we chose for the 
initial magnetic and bulk velocity fluctuations, ac and 
ctr are both statistically very close to zero, even though 
they are not actually zero anywhere in the (t, y) plane. 

A non-zero resistivity has been introduced in order to 
guarantee a satisfactory conservation of the total energy, 
with no claim to model any realistic physical process. Its 
value has been empirically fine-tuned by running differ¬ 
ent simulations (see TableJ^. Eurther details about this 
point will be provided in Sec. 


3. RESULTS 

We performed nine different simulations. Their main 
parameters are listed in Table A label is assigned to 
each run in the first column, while in the other columns 
we report, from left to right: the spatial resolution. Ax 
(= A^), the length of the simulation box, Lbox, the value 
of the resistivity coefficient, 77 , and the number of ppc. 

















4 


FRANCI ET AL. 


Run A employs the best spatial resolution, 
Ax = = 0.125 (ip, and the highest number of 

particles, i.e., 8000 ppc, corresponding to more than 
3 X 10^^ particles in the whole simulation domain. The 
resistive coefficient has been fine-tuned and set to the 
value 77 = 5 X 10“^, in units of A7tuj~^. In Subectionf 


action |3.l| 
tive anal- 


and |3.2| we will provide a detailed and quantitative anal 
ysis of the data produced by this run. The remaining 
simulations were performed in order to validate these 
results and to investigate the effects of the number of 
ppc (Runs B-E, see Table Q, the resistivity (Runs E-G), 
and the spatial resolution (Runs H-I). Their results will 
be discussed later, in Section 

3.1. Temporal and spatial evolution 

In Eig. the time evolution of a few quantities is 
shown up to 500 The initial non-linear time associ¬ 
ated to the maximum injection scale, i.e., ^ 0.2 

can be estimated as ^nl ^ ^ 2011 “^, and 

corresponds to the minor ticks of the x axis. The total 
length of the simulation is approximately 25 ^nl- 

In the first panel, from top to bottom, we report the 
rms out-of-plane current density, Jy, and the rms out-of- 
plane vorticity, a;||. The current increases quite rapidly, 
attains its maximum value just before t r\j 200 ftp ^ and 
then slowly decreases. Since it represe nts a good indi¬ 


cator of th e level of turbulent activity (Mininni & Pou- 
quet 2009), we choose to perform a detailed analysis at 


t ^ 200 ftp ^ , when the turbulence is expected to be fully 
developecl. A vertical black dotted line marks this time in 
all four panels. The vorticity also increases quite rapidly, 
reaching an earlier and lower maximum value, and then 
it decreases extremely slowly. 

The second panel of Eig. [l] shows the rms perpendicu¬ 
lar and parallel magnetic fluctuations (red lines), 
and the rms perpendicular uj_ and parallel u\\ velocity 
fluctuations (blue lines). Perpendicular and parallel com¬ 
ponents are drawn with solid and dashed lines respec¬ 
tively. B± exhibits a small increase until t ^ 40 0“^, and 
then it decreases quite smoothly. On the other hand, Uj_ 
decreases with a similar trend, but without showing any 
initial growth. This indicates that the turbulence is fed 
by the perpendicular components of both the magnetic 
and the velocity fluctuations, whose energy decreases 
slowly and sustains the cascade for the whole evolution. 
Contextually, the parallel components of both the mag¬ 
netic and the velocity fluctuations rapidly originate from 
compressive effects, remaining much smaller than their 
perpendicular counterparts throughout the simulation. 

In the third panel of Eig. we report the space- 
averaged parallel and perpendicular proton tempera¬ 
tures, normalized to the initial value Tq, {Tp\\/To) and 
{Tp±/To) respectively, together with the space-averaged 
proton temperature anisotropy, (Ap) = {Tp_i/Tp\\). We 
recall here that Ty and T±_ are (defined with respect to 
the local magnetic field. Since {Ap) = 1 is imposed at 
t = 0 , Tpii and Tp± share the same initial value, Tq. In 
the very first part of the evolution, the former shows 
a little and sudden decrease, after which both increase 
with almost the same rate. The parallel and perpendic¬ 
ular energy gains, at the end of the simulation, i.e., at 

respectively. 


t = 500 ft^ are approximately 6% and 



0 100 200 300 400 500 

t/Qp 


Fig. 1. — Time evolution of space-averaged quantities. From top 
to bottom: the rms out-of-plane current density, J y, and the rms 
out-of-plane vorticity, a;|| {first panel); the rms perpendicular B± 
and parallel Sy magnetic fluctuations, and the rms perpendicu¬ 
lar u± and parallel ity velocity fluctuations {second panel); the 
mean values of the normalized perpendicular and parallel proton 
temperatures, Tp_L/To and T^y/To, and of the proton temperature 
anisotropy, Ap {third panel); the mean values of the normalized 
cross helicity, ac , and of the normalized residual energy, {fourth 
panel). In all panels, a vertical black dotted line marks the time 
of the maximum turbulent activity, i.e., t = 20012^ . 

This small excess of perpendicular energy quickly arises 
within ^ 2 ^nl , and it is preserved throughout the simu¬ 
lation, with the temperature anisotropy reaching a value 
{Ap) 1.02 in correspondence of the maximum turbu¬ 
lent activity, and then remaining quite constant until the 
end of the simulation. A detailed discussion about the 
proton heating will be further provided in Subsection |4.2[ 
Lastly, in the bottom panel of Eig. we show the 
space-averaged values of the normalized cross helicity, 
ac , and of the residual energy, cfr^ (see Eq. © and 0 ). 
The former is very close to zero at the beginning ofme 
simulation (as a result of the initially imposed random 
phases spectra), and it tends to maintain this value un¬ 
til the end. The latter, instead, decreases from zero to 
about —0.3 in very few non-linear times ^nl, showing a 
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Fig. 2.— Contour plots of six different quantities on the {x^y) plane at t = 200the amplitude of the perpendicular magnetic 
fluctuations, {upper-left panel), and of the perpendicular velocity fluctuations, {upper-right panel), the out-of-plane current density, 
J|| {middle-left panel), and vorticity, a;|| {middle-right panel), the proton temperature variation normalized to the initial temperature, 
ATp/To {bottom-left panel), and the proton temperature anisotropy, Ap {bottom-right panel). 
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Fig. 3. — Top panel: PDFs of the increments of the perpendicular 
magnetic field component By along x at t = 200 correspond¬ 

ing to 27Tdpli = 0.4, 27Tdpli = 2 and 2ixdplt = 8 (from top to 
bottom). In each panel, a Gaussian function with the same vari¬ 
ance is plotted with a dashed line as a reference. Bottom panel: 
Excess kurtosis of the same quantity, computed at the same time. 


global excess of the magnetic energy over the kinetic en¬ 
ergy. These asymptotic values are reached very quickly, 
as a consequence of the relaxation from the initial ran¬ 
dom relative orientation of the velocity and the magnetic 
fluctuations towards a strongly aligned state. Despite 
the steady time evolution of their space-averaged val¬ 
ues, both Gc and gr appear very patchy when looking 
at the spatial distribution throughout the 2D computa¬ 
tional domain (not shown), exhibiting quite a wide excur¬ 
sion from —1 to 1 between different albeit close regions. 

Summarizing the time evolution of all the above- 
mentioned quantities, we can divide the evolution of the 
system in three different stages: 

1. a rapid re-adjustment and relaxation of the initial 
conditions, occurring within t < 40 ^ 2tNL 


2. the onset of a turbulent cascade, fed by the perpen¬ 

dicular magnetic and velocity fluctuations, involv¬ 
ing larger and larger scales on times of the order of 
t ^ 200 ^ 10 

3. a deacaying phase with slow and smooth variations 
of all rms quantities, during which the turbulence is 
fully developed and further sustained until at least 
t ^ 50011“^, corresponding to ^ 25 ^nl- 

Fig. |2] shows isocontours of six different quantities 
in the whole simulation domain, all computed at t = 
200 In the upper-left panel, we report the magni¬ 
tude of the perpendicular magnetic fluctuations, 
showing the presence of coherent structures in the mag¬ 
netic field, i.e., vortices and magnetic islands, embedded 
in a much more chaotic environment where stretched 
and twisted shapes emerge. In the upper-right panel, 
the magnitude of the perpendicular velocity fluctuations, 
is shown to exhibit qualitatively the same kind 
of structures, but with lower intensity and much lower 
gradients. In some regions, high values of corre¬ 
spond to high values of while in other regions the 

opposite situation holds. In the middle-left panel, we 
show the out-of-plane current density, Jy = (V x -B)||. 
Many thin current sheets form, since the very first phase 
of the evolution, mostly around and in-between vortices. 
Once formed, each current sheet is quickly disrupted into 
smaller and smaller pieces, contributing to the genera¬ 
tion of smaller-scale structures. At the time of maximum 
turbulent activity, this results in the articulated pattern 
shown here. In the middle-right panel, the out-of-plane 
vorticity, a;|| = (V x n)||, is shown. It exhibits a struc¬ 
ture similar to Jy, with many thin layers, whose shape 
is however much more defined and clean in respect to 
Jy. Peaks of cjy and peaks of Jy occupy approximately 
the same regions. In the bottom left panel, we report 
the proton temperature variation in respect to the ini¬ 
tial proton temperature, AT^/Tq = {Tp — To)/To, where 
Tp = {2Tp^ + Tpy)/3 is the average proton temperature 
measured at t = 200 Regions where AT^ is locally 
both negative or positive are clearly present, and a re¬ 
sulting global proton temperature enhancement can be 
observed, as already inferred from Fig. Interestingly, 
areas where a proton temperature enhancement o ccurs 
are located in th e vicinity of current sheets (cf., (|Ser- 
vidio et ^|20I2 )). A more detailed analysis shows that 
strong currents exhibit a complex evolution, which in¬ 
volves splitting/dissociation and leads to a relevant pro¬ 
ton energization. 

In the bottom right panel, the proton temperature 
anisotropy, Ap, is shown. We observe a wide excursion 
between very close areas, the perpendicular proton tem¬ 
perature Tp^ ranging from about half and almost twice 
the parallel one. Therefore, there is a strong local reshap¬ 
ing of particle distributions, l eading to both perpe ndicu- 
lar and parallel anisotropies ( Seryidio et al.|[2M4|). Nev¬ 
ertheless, as inferable from Fig. the relative difference 
between {Tp±) and {Tp\\) is about 2% at t = 200 
meaning that globally no preferential enhancement along 
the perpendicular or parallel direction is achieved. 

The small-scale coherent structures which have 
emerged by the time of maximum turbulent activity. 













Kinetic plasma turbulence at proton scales 


7 



k^dp 


Fig. 4.— Power spectra of the perpendicular magnetic and ve¬ 
locity fluctuations, (red solid line) and (blue solid line), 
respectively. Power laws with different spectral indices are addi¬ 
tionally shown in black dashed lines as a reference. 


already observed in Fig. can be related to the phe¬ 
nomenon of intermittency, since they are able to induce 
departures from self-similarity and enhanced dissipation. 
In order to look for intermittency in our data, we exam¬ 
ine the non-Gaussian behavior of the probability density 
function (PDF) of a MHD primitive variable. In particu¬ 
lar, we compute the PDFs at t = 200 by taking incre¬ 
ments of one of the perpendicular magnetic field compo¬ 
nents, i.e., By^ along the other perpendicular direction, 
i.e., X, for three different spatial separations, In the 
three top panels of Fig. we show three PDFs, com¬ 
puted for = 0.4, ^ich is approximatively in the 

middle of the inertial range (upper panel), 2'Kdpji = 2, 
which is the scale corresponding to the ion spectral break 
(middle panel) and 27rdpl£ = 8, which is well inside the 
kinetic range (bottom panel). A Gaussian function with 
the same variance is plotted with a dashed line in each 
panel as a reference. The distribution of magnetic fluctu¬ 
ations is clearly different at different scales: it is closer to 
a normal distribution at very large scales, it shows a sig¬ 
nificant deviation at intermediate scales, and it displays 
very extended tails at small scales. In order to quan¬ 
tify the level of intermittency, we compute the fourth 
central moment K (or kurtosis) of the distributions. In 
the bottom panel of Fig. we show the excess kurto¬ 
sis A"' = — 3, computed from the increments of By 

along X, as a function of 2'Kdpji^ again at t = 200^4“^. 
This quantity is clearly very close to zero up to the injec¬ 
tion scale, i.e., k^dp < 0.2, and the it steadily increases 
through the inertial range and down to sub-proton scales. 


3.2. Spectral properties 

Since the small-scale structures shown in Fig. [^exhibit 
random orientations, and therefore the two-dimensional 
spectra of all fluctuations can be assumed statistically 
isotropic, we can perform a quantitative analysis of 
the turbulent cascade by computing the omnidirectional 
spectra. These are defined as 

P^{ki)=5'i>\ki)/k^= Y. (5) 

\h^\=k^ 



0.1 1.0 


/cy dp 

Fig. 5. — Top panel: Power spectra of the perpendicular mag- 

5/3 

netic ffactuations, ^_l, compensated by . Each of them is the 

time average in the interval [t — 20r2p^, t + 20r2p^], where i is 
the time reported in the legend. Note that all of them have been 
suitably rescaled for the sake of clarity, so that they did not over¬ 
lapped each other. Horizontal dotted black lines are additionally 
shown. Bottom panel: The same as in the upper panel but for the 

perpendicular velocity ffuctuations, u±, compensated by k^‘^. 


where T are the Fourier coefficients of a given quantity 
T, and 6'^{k±) is the amplitude of the fluctuation T at 
the scale k_\_. 

In Fig. we show the spectra of the perpendicular 
magnetic and velocity fluctuations, drawn with red and 
blue solid lines respectively, at t = 200^4“^. We clearly 


observe two power-law ranges, separated by a smooth 
spectral break at a scale of th e order of the the proton 
inertial length, k± dp 2. In jLetter 1| we showed the 
spectra of the total magnetic and velocity fluctuations, 
which exhibit a very similar behavior, since the perpen¬ 
dicular components are the dominant ones for both field. 

In particular, in the inertial range the spectrum of 
the perpendicular magnetic fluctuations follows a Kol¬ 
mogorov power-law scaling over a full decade in 

wavenumber, approximately between k±dp = 0.1 and 
2. Simultaneously, the perpendicular proton bulk ve¬ 
locity fluctuations exhibit a less steep slope, with an 

Iroshnikov-Kraichnan scaling of kj_ ' , over a little less 
than a decade, in the range 0.2 < k±dp < 1. More¬ 
over, an excess of magnetic energy over kinetic energy is 
observed, coherently with the negative value of the nor¬ 
malized residual energy ctr already shown in the bottom 
panel of Fig. 

At kinetic scales, the spectra of both fields steepen, 
due to the presence of both kinetic and dissipative (re¬ 
sistive) effects. The spectrum of u± quickly drops with 
an exponential trend above k±dp ~ 1, until it clearly sat¬ 
urates to the noise level due to the flnte number of ppc. 
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Fig. 6. — Spectra of the total current density, J (grey line), of 
the density fluctuations, n (purple line), of the magnitude of the 
magnetic field, \B\ (orange line), and of its parallel component, Sy 
(red dot-dashed line). Power laws with different spectral indices are 
additionally drawn in black dashed lines, as a reference. 

corresponding to the spectrum at t = 0. The spectrum 
of the magnetic fluctuations, on the contrary, continue 
to follow a power-law scaling also at sub-proton scales, 
although with a steeper spectral index, of the order of 
—3. For k_\_dp > 10, Pb^ does not show an exponential 
damping, as one would expect for resistive dissipation, 
but a small increase instead, since the adopted resistive 
coefficient is slightly smaller than the optimal value. 

As discussed, the maximum level of turbulent activity 
occurs at t ^ 200 , which is about ten times the initial 

nominal nonlinear time ^nl- This can be explained with 
the fact that at t = 0 we inject energy through several 
modes within the range [/cq, where ko is the largest 
scale corresponding to the computational box size, i.e., 
ko = 27r/(256dp) ^ 0.025The nominal nonlinear 
time tNL|t=o is different for each mode, being longer for 
lower /c-vectors. As the system evolves, the injection 
scale gets larger and larger and most of the initial modes 
are involved in the development of the turbulent cascade 
at t = 200^2“^. Since modes with lower ks keep feeding 
energy at large scales even afterwards, we expect turbu¬ 
lence to be still sustained also at later times. 

In Fig. we show the spectra of the perpendicular mag- 

netic fluctuations, compensated by k_^ (top panel), and 
the spectra of the perpendicular velocity fluctuations, 

compensated by k^/‘^ (bottom panel), computed at reg¬ 
ular intervals of 50^2“^, from the maximum of turbulent 
activity to almost the end of the simulation. Here, the 
spectrum at a given time i is the time-average between 
five different spectra corresponding to t, t ± 10 ftp^ and 
t ± 20 ^p^. The power-law scaling for both the magnetic 
and the velocity fluctuations are very well maintained, 
over about the same range, at all times t > 200 in¬ 
dicating that the turbulence decays in a self-similar way. 
Note that spectra corresponding to different times have 
been slightly shifted along the vertical axis, in order to 
avoid overlapping. 

In Fig. we show the spectra of the magnitude of 
the magnetic field, \B\ (orange), of its parallel compo- 



k±dp 

Fig. 7. — Spectrum of the perpendicular electric field, Ej_ (green 
solid line) and energy associated to the different terms of Eq. 

(the term containing the resistive coefficient is negligible). A power 
law with a spectral index of —0.8 is also drawn with a dashed black 
line, as a reference. The shaded grey region marks the range where 
numerical effects strongly affect the shape of Pej^ ■ 

nent, (red dot-dashed), of the density fluctuations, 
n (purple) and of the total current density, J (grey). 
The density and the parallel magnetic fluctuations are 
strongly coupled beyond k^dp ^ 2. In the MHD range, 
they exhibit a flat spectrum, which is approximately an 
order of magnitude smaller than the one of the perpen¬ 
dicular magnetic fluctuations. Therefore, the large-scale 
activity has little contribution from compressible fluctu¬ 
ations - although they can still play a significant role in 
the dynamics of the out-of-plane components - and the 
magnetic compressibility, i.e., the ratio of parallel to to¬ 
tal magnetic fluctuations, is also negligible at small /cs. 
Both spectra steepen at sub-proton scales, following a 
clean power-law scaling with a spectral index of —2.8. 
Note that their relative power level with respect to other 
fields’ spectra increases, with Pb|| (and also P\b\) becom¬ 
ing comparable with the ^ectrum of the perpendicular 
component, Bj^ (cf. Fig. ffl. 

The spectral shape of the current density, J, can be 
understood by recalling that J = V x S. A simple order- 
of-magnitude estimate of its perpendicular and parallel 
components gives J± oc k±B\\ and Jy oc kj_Bj_^ respec¬ 
tively. Therefore, in the inertial range, where the mag¬ 
netic activity is dominated by the perpendicular fluctu¬ 
ations, the spectrum of J is determined by its parallel 
component J\\ and this results in the observed spectral 

index of +1/3, since oc with Pb^ cx: k^^^ . 

On the contrary, as already discussed, Pb,, and Pe^ are 
comparable at sub-proton scales, both showing a power- 
law scaling, with spectral indices of —2.8 and —3, respec¬ 
tively. The corresponding components, Jy and J^, are 
of the same order and exhibit a similar scaling, therefore 
Pj Pj\\ corresponding scaling is in- 

between oc k~^-^ and oc k~^. The change in the spectral 
slope of Pj at kj_dp 2 provides a further evidence of a 
spectral break in the magnetic field spectrum at proton 
scales. 

Finally, the spectrum of the perpendicular electric fluc¬ 
tuations is reported in Fig.[^ as a green line. We choose 
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Fig. 8. — Spectrum of the perpendicular electric field, as in 
Fig. 1^ for the case with = 0 (Vpe = 0). 


to pay particular attention to the electric field for three 
main reasons. Uppermost, it is expected to exhibit the 
most complex spectrum, since it contains the contribu¬ 
tions of four terms having different relative importance in 
different ranges of scales. Secondly, it is the quantity that 
is mostly affected by both numerical effects and particle 
properties, so its behavior needs to be analyzed carefully, 
especially at small scales. Lastly, no consistent observa¬ 
tional data about the properties of E are available yet, 
so making predictions about the shape of its spectrum 
can be relevant for future analysis. We recall here that, 
starting from the Vlasov-fluid equations and assuming 
that the electrons act as a massless, charge-neutralizing 
fluid, the electric field can be computed from the gener¬ 
alized Ohm’s law as 



(6) 


^MHD 


In Fig. together with Pe^ obtained from numerical 
data, we also report the energy associated to the first 
three terms of Eq. (§, computed a-posteriori and drawn 
with cyan, magenta and black dot-dashed lines respec¬ 
tively (the contribution from the resistive term is neg¬ 
ligible at all scales, since the resistive coefficient 77 is 
5 X 10“^.) At large scales, Pe^ is clearly dominated 
by the MHD term, F^mhd, which is essentially perpen¬ 
dicular to ^ 0 , since its leading contribution comes from 
u± X Bq {uj_ X B\\ and x B± are both of the second- 
order in the fluctuations). Therefore, Pe^ follows strictly 
the spectrum of the perpendicular velocity fluctuations 
(cf. Fig. [^. Approximately at k±dp ^ 0.5, these two 
spectra decouple, since the second and third terms of 
Eq. ^ start contributing. 

We can accurately analyze the Hall term, E^Hain by 
considering its perpendicular and parallel components 
separately. The former is of the first-order in the fluctu¬ 
ations, being led by {k± x S||) x Bq (other contributions 
are quadratic in -By and B±, and therefore negligible). 
On the contrary, the latter is only of the second-order 
in the fluctuations, having the only contribution from 
{k±xB\\) X B±. Therefore, we expect the Hall term to 


be negligible at large scales, where J± is small, and to 
exhibit a power-law behavior at small scales, with spec¬ 
tral index ^ — 0.8 following from Enaii cx: qPB„. In- 
deed, this is what we observe in Fig. (compare the 
magenta dot-dashed line with the reference dashed black 
line). The electron pressure gradient term, Bpe, has only 
perpendicular components by construction (our 2D com¬ 
putational domain is perpendicular to Bq). In the iner¬ 
tial range, it is of course negligible, the spectrum of the 
dens^ fluctuations being essentially flat (compare with 
Fig.[^. On the contrary, at small scales, it is expected to 
give a contribution Pvn oc k\Pn^ which has exactly the 
same slope as the contribution from the Hall term, since 
the spectra of the density fluctuations and of the paral¬ 
lel magnetic fluctuations have the same spectral index of 
—2.8 at sub-proton scales. This is indeed what we ob¬ 
serve in Fig. where the contribution from the electron 
pressure gradient term is drawn with a black dot-dashed 
line. We would expect a similar behavior for Pe^ at sub¬ 
proton scales, i.e., a power law with a spectral index of 
^ — 0 . 8 , and we observe a hint of a similar scaling in the 
range 2 <k^dp <1. 

The spectrum of the electric field fluctuations is the 
most affected by numerical effects among all the consid¬ 
ered spectra, since the computation of E involves both 
other fields {u and B) and derivatives (V x B and Vn), 
as shown by Eq. We already noticed that Pb suffers 
from an accumulation of energy at small scales, which is 
only a numerical artifact, and so does Pj (cf. Fig. [^. 
Moreover, derivatives in the numerical code are com¬ 
puted as finite differences, thus they are not able to re¬ 
cover very precise quantities at very small scales. For all 
these reasons, we cannot extract any robust information 
about the spectrum of the electric field at high wavevec- 
tors. In order to emphasize this, we choose to mark the 
“non-safety area”, which we estimate as k±dp > 7, with 
a shaded gray region in Fig. 

Under particular conditions, i.e., = 0, one can ob¬ 

tain a better defined scaling for the electric field. In this 
case, the electron pressure CTadient term, Bpe, in the 
Ohm’s law is zero (see Eq. Q). Consequently, the level 
of the electric field spectrum at small scales is higher, 
since it is supported only by the Hall term Bnaii- This 
can be seen in Fig. where we show the same analysis 
of the electric field spectrum as in Fig. but for a case 
with l3e = 0 and all the other parameters set as in Run 
A. The electric field spectrum is now the sum of only 
two main contributions, Bmhd and BHaii- No qualita¬ 
tive changes are introduced with respect to the case with 
a finite electron temperature; the former term dominates 
the spectrum at MHD scales, while the latter is respon¬ 
sible for the flattening of the electric field at ion scales. 
The important difference with respect to Run A, is that 
Pej^ displays now a well defined power-law slope with 
an index of —0.8, consistent with the expectation. We 
expect the same slope also for the finite case of Run 
A, in the absence of the numerical limitations discussed 
above. 

4. ROLE OF NUMERICAL PARAMETERS 
4.1. Spectral properties 

As mentioned, a small numerical resistivity, 77 , has been 
implemented in all runs (see Table [^. A proper level of 
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dp 

Fig. 9. — Power spectra of the perpendicular magnetic fluctua¬ 
tions, S_L, for simulations with different values of the resistivity, 
i.e.. Run A, (ry = 5 x 10“"^, black line). Run F (77 = 1 x 10“"^, red 
line) and Run G (77 = 1 X 10“^, blue line). 

resistivity is mandatory in order to prevent an accumu¬ 
lation of the energy in magnetic fluctuations at small 
scales, due to numerical errors. Runs F and G have been 
used to fine-tune the resistivity coefficient, starting from 
an order-of-magnitude estimate, and are characterized 
by the values 77 = 1 x 10“^ and 1 x 10“^, respectively. In 
Fig.f^ we show the corresponding spectra of the perpen- 
dicu& magnetic fluctuations, at t = 200 in 

comparison with Run A {rj = 5 x 10“^). For the setting 
adopted, the dissipative scale for the under-resistive case 
(Run F) can be estimated as k^isdp ^ 35, i.e., smaller 
than the scale corresponding to the employed resolution. 
As a consequence, this value of the resistivity is not high 
enough to remove the energy excess at small-scales, as 
also demonstrated by the shape of the spectrum of Run 
F in the sub-proton range. The over-resistive simula¬ 
tion (Run G) corresponds to the opposite case, where 
^dis dp ^ 6 , thus well inside the range of wavevectors re¬ 
solved in the simulation, making Pbj_ decreasing expo¬ 
nentially below the ion-break, a clear indication of a too 
strong dissipative damping at sub-proton scales. Lastly, 
the intermediate case, i.e.. Run A, is expected to intro¬ 
duce a dissipation scale k^is dp ^ 10 , allowing for the best 
description of the spectrum of B±. Indeed, this is ob¬ 
served to follow a power-law scaling with a spectral index 
of —3 for roughly a decade after the break (cf. Fig. [^, 
in good agreement with solar wind spectra from obser¬ 
vational data. Therefore, we decided to adopt 5 x 10“^ 
as the optimal value for the resistive coefficient 77 . 

As a further confirmation for the adequacy of our 
choice for 77 , in the insert of Fig. we also show the time 
evolution of the total energy f, normalized to its initial 
value, for the same three simulations. In all three cases, 
the total energy stays constant for t < 70(1“^, while a 
different behavior is observed at later times. When the 
resistivity is too low (RunF, red line) £ grows signifi¬ 
cantly, due to the inefficient control of energy at small 
scales, and such an increase is already of the order of 
^ 1 % at t = 200 On the contrary, when 77 is too 

high (Run G, blue line) the action of resistivity is too 
strong. Note that part of the energy subtracted from the 



k± dp 


Fig. 10.— Ratio between the omnidirectional spectra of per¬ 
pendicular electric ffuctuations, Pej^ , and perpendicular magnetic 
ffuctuations, Pbj^ 5 for simulations with different resolution: Run 
A, with Ax = 0.125 dp (black line). Run H, with Ax = 0.25 dp 
(blue line), and Run 1, with Ax = 0.5dp (red line). 

magnetic fluctuations would also go into electron heat¬ 
ing but, since the hybrid approximation does not provide 
an evolution for the electron temperature, such energy is 
not taken into account, and is then lost by the system, 
resulting in a net decrease of £. The value 77 = 5 x 10 “^ 
is the one which better ensures the conservation of the 
total energy, with a relative difference between the begin¬ 
ning and the end of the simulation, i.e., t = 50011“^, of 
about 0.3%. Also note that, although the shape of Pbj_ 
at sub-proton scales is quite strongly affected by the re¬ 
sistivity, the power-law scaling in the inertial range and 
the position of the spectral break are not, assuring the 
reliability of the spectra shown in Fig. Hand Fig. 

We also investigated the stability of omnidirectional 
spectra versus the spatial resolution. This was done by 
varying Ax, while keeping fixed the number of grid points 
(the total box length is then larger for larger grid spac¬ 
ing) , as well as the amplitude of the initial magnetic fluc¬ 
tuations, Run H and I implement Ax = 0.250 

and 0.500dp, respectively (see Table [T|) . For both these 
runs, the value of the resistivity coefficient was suitably 
rescaled in order to get dissipation at the proper scales 
(note that 77 oc Ax under these conditions). 

In Fig.[^ the ratio between the omnidirectional spec¬ 
tra of the perpendicular electric and magnetic fluctua¬ 
tions, Pej_IPbj_^ is compared between Run A, H and I. 
For all the three runs, this ratio exhibits the same scal¬ 
ing in the inertial range, following a power law with a 
spectral index of 1/6. This is a direct consequence of 
the different scaling for the magnetic held (—5/3) and 
the velocity (—3/2) in the ideal MHD regime - where 
also Pe^ ^ Puj_ - leading then to a spectral index 
of —3/2 -j- 5/3 = 1/6 for their ratio. As the other 
terms in the generalized Ohm’s law became important 
at ion scales, the Pe^IPb^ ratio increases significantly 
at smaller ks. Interestingly, increasing Ax from 0.125 
(blue) to 0.250 dp (black) does not produce a change in 
the scale at which Pe^ exceeds Pbj_^ and this break is 
observed to occur at k±dp 2 in both cases. Moreover, 
the two curves exhibit similar slopes in the sub-proton 
ranges. This is a confirmation that the estimate of 77 
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Fig. 11.— Zoom of the power spectra of the perpendicular mag¬ 
netic, electric and proton bulk velocity fluctuations (drawn with 
red, green and blue lines, respectively) at small scales, where the 
contribution of numerical noise is not negligible. Lines with differ¬ 
ent shades of the same color correspond to simulations with differ¬ 
ent amounts of ppc, ranging from 500 to 8000, with darker colors 
being associated with a higher number of particles. 


for the two simulations was correct and that the raise 
in Pej_IPbj_ is physical and well captured by the runs. 
On the other hand, when employing a lower resolution, 
Ax = 0.500 (red), the break seems to occur at slightly 
larger scales. This is likely a consequence of the reduc¬ 
tion of the resolution at small scales: in Run I, the break 
and the dissipative scale are not well separated in Fourier 
space, so that subtracting energy at the smallest scales 
via resistive dissipation also affects the shape of the spec¬ 
tra around k_\_ dp 2, where the break occurs. This is an 
evidence that the scaling for the spectra discussed and 
shown in Fig. continue to hold at lower spatial reso¬ 
lution, but also that Ax > 0.500 dp is not sufficient to 
properly explore the physical behavior at sub-ion scales. 


Finally, the importance of employing a high number 
of particles was investigated, by keeping all the same 
parameters as for Run A except for varying the num¬ 
ber of ppc from 4000 to 500 in steps of a factor of 2 
(Runs from B to E in Table [^. In Fig. 11, the spec¬ 
tra of the perpendicular velocity, magnetic and electric 
fluctuations are reported with lines in different shades of 
blue, red and green respectively, corresponding to simu¬ 
lations with different numbers of ppc ranging from 500 
(lighter color. Run E) to 8000 (darker color. Run A). In¬ 
creasing the number of pec from 500 (^ 2.0 x 10^ total 
particles) to 8000 (^ 3.3 x 10^^ total particles) results in 
an decrease of the noise at small scales of more than one 
order of magnitude for the spectrum of the perpendicu¬ 
lar velocity fluctuations. On the other hand, the trend of 
the spectrum up to the proton inertial length and slightly 
above is not affected, and the curve corresponding to our 
most accurate simulation (Run A) overlaps with the one 
with 4000 ppc (Run B) up to k±dp ^ 4. On the contrary, 
the spectrum of the perpendicular magnetic fluctuations 
is only barely affected by the numerical noise when the 
number of particles is sufficiently high, since the curves 
for 4000 and 8000 ppc are almost indistinguishable, prov¬ 
ing that the number of ppc employed in Run A is suf¬ 
ficient to get reliable results for Pb^ up to k_\_dp ~ 10 . 


Lastly, the spectrum of the perpendicular electric fluc¬ 
tuations shows a dependence on the number of particles 
at scales k_\_dp > 4, but the curves for 4000 and 8000 
ppc are quite close to each other even at smaller scales. 
However, as mentioned, Pe^ is influenced by different 
sources of numerical noise, and all contribute in affect¬ 
ing the spectrum at small scales. We would like to stress 
that the evaluation of the noise due to the finite number 
of ppc only represents a lower limit of the overall noise, 
and therefore our previous choice of marking a shaded 
grey area for k±dp> 7 in Fig. 0 is not in contrast with 
these results. 


4.2. Proton heating 

Eig. shows that some particle heating is observed 
during the turbulent activity. Some care must be used in 
the interpretations of this result, since it may be signifi¬ 
cantly affected by some of the numerical settings. There¬ 
fore, we carefully consider the properties of the proton 
heating in this subsection. 

Resistivity is observed to play a fundamental role 
in determining the overall proton heating, ATy ^ = 
(^||,±) ~ ^ 0 , and the proton temperature anisotropy, Ap. 
In Fig.[^ we show the time evolution of the perpendicu¬ 
lar Tp±_ and the parallel Tp\\ proton temperature, in solid 
and dashed lines, respectively, corresponding to different 
values of the resistive coefficient rj (Run A, E and G of 
Table JTl). The time evolution of T^n is observed to be 
not affected almost at all by the resistivity, showing an 
early decrease up to t ^ 50f)“^ and then an increase 
with an almost constant rate, as was already shown for 
Run A in Fig. The situation is different for Tp±^ since 
its behavior for different values of 77 is the same only 
during the initial readjustment of the system, while it 
starts to differ after t ^ 40^4“^. At later times, Tp±_ ex¬ 
hibits a growth rate very similar to that of T^n for Run A 
(black), so no preferential perpendicular or parallel heat¬ 
ing is observed. In particular, AT^^/To at t = 200 ftp^ is 
about 3.5%, while the corresponding ATpy/To is about 
2%. When 77 is lower (Run F, red), Tp±_ grows with a 
much faster rate than '^pW for t > 50 ^, generating 

a strong preferential heating in the perpendicular direc¬ 
tion, Tpj_ being about 8 % greater than the initial value 
at t = 200^4“^. On the contrary, when 77 is higher (Run 
G, black lines), Tpj_ grows much slower, being overcome 
by Tp|| just before t = 200 0 “^, leading to < 1 at later 
times. The amount of perpendicular heating observed is 
then significantly related to the presence of an excess of 
fluctuations at small scales, and can be then therefore 
largely overestimated, or underestimated, if an incorrect 
value of the resistivity is adopted. 

The perpendicular proton heating is also found to be 
strongly affected by the number of particles employed. 
In Fig.|T^ we report the ratio ATp^/To, where ATp± = 
{Ppi.) — To is computed at t = 200 versus the num¬ 
ber of ppc for Runs from A to E (see Table 0. At this 
time, ATpj_/To is clearly higher when employing a lower 
number of ppc. In particular, while a good convergence 
towards a value of ^ 3.4% is observed when increas¬ 
ing the number of ppc from 4000 to 8000, this quantity 
clearly diverges when few particles are employed, reach¬ 
ing 4.8% for 500 ppc. Moreover, the difference in Tp± 
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Fig. 12. — Time evolution of the parallel and perpendicular pro¬ 
ton temperature, T^n and Tpj_, respectively, normalized to the ini¬ 
tial common value, Tq. The evolution is here shown for different 
values of the resistive coefficient (see Table [^. Solid and dashed 
lines are used for and Tp||, respectively. 



Fig. 13. — Perpendicular proton heating, ATpj_/To, computed 
at t = 200 versus the number of ppc employed, ranging from 
8000 (Run A) to 500 (Run E). 


between Run A and Run E tends to further increase at 
later times. The main result of this analysis is that the 
use of a high number of ppc is mandatory when trying to 
give an estimate of Tp±, which could be largely overes¬ 
timated otherwise. On the contrary, the parallel proton 
temperature is found to be largely independent from the 
number of particles (the relative difference between Run 
A and Run E is lower than 0.1% at t = 20011“^). 

To conclude our analysis, we find that the spatial res¬ 
olution Ax does not seem to affect significantly the pro¬ 
ton heating, provided that the value of the resistivity is 
suitably set as discussed in Subsection |4.1[ Differences 
between Tp±_ at t = 200 for Runs A, H, and I are 
less than 0.6%. The dependency of T^n on the spatial 
resolution is also negligible. 

5. SUMMARY AND CONCLUSIONS 

In this work, we have presented properties of turbu¬ 
lence in a magnetized collisionless plasma by means of 


two-dime nsional hy brid PIC simulations, extending the 
results of |Letter~T Remarkably, our simulations imple¬ 
ment a high number of collocation points (2048 x 2048) 
and a very high number of particles (up to 8000 ppc), 
covering a large simulation domain (not less than Lbox = 
256 dp) with a fine spatial resolution. This enables to 
self-consistently describe the evolution of turbulence over 
three orders of magnitude in wavevectors, and to fully 
capture its transition from fluid-hke MHD scales to ki- 
netic sub -ion scales, by using a single simulation (see 
Letter I). 

T'he adopted initial conditions consist of balanced and 
equipartitioned magnetic and velocity fluctuations, i.e., 
with zero cross helicity and zero residual energy. The 
onset of a turbulent cascade appears quite early during 
the simulations (t ^ 200 corresponding to approxi- 
matively 10 nonlinear times ^nl), h e., when most of the 
initial modes have started to partake into the cascade. 
In physical space, the activity of turbulence is character¬ 
ized by magnetic field coherent structures, vortices, and 
strong and localized current sheets at smaller scales. 

Generation of coherent structures associated to inter- 
mittency is observed as turbulence evolves through MHD 
to sub-proton scales. PDEs of increments of a per¬ 
pendicular component of the magnetic fluctuations at 
t = 200 exhibit a deviation from the normal dis¬ 
tribution at all scales. This is small in the inertial range, 
becoming larger in correspondence of the spectral break 
at k± dp 2, while at k_\_ dp 2 the PDE has a much 
leaner shape with long non-zero tails. The corresponding 
excess kurtosis confirms this behavior. It is very small 
around the injection scale, since part of the MHD range 
fluctuations still acts as an energy reservoir for turbu¬ 
lence at t = 200 while it increases through the in¬ 
ertial range. Moreover, we observe a further increase at 
smaller scales. Observational data give no firm results 


(1 Alexandrova et al.||2008a Kiyani et al.| 

|2009| 

, 2013 

Wu 

et al.|2UI3 Chen et al.|2(Jl4|). JNevertheless, wl 

len shown. 


tosis at smaller scales (e.g., 

Dmitruk & Matthaeus 

||2QQ6l 

Wan et al.|20I2 

Karimabad 

i et al.|2U13 Wu et al.l* 


ties, two clear distinct turbulent regimes are observed. 
At larger scales, the magnetic field follows a Kolmogorov 
—5/3 power law, while the velocity has a spectral index of 
—3/2, which is characteristic of a Iroshnikov-Kraichnan 
turbulence. An excess of magnetic energy with respect 
to the kinetic energy is observed throughout the inertial 
range. The two different scalings for the magnetic and 


(iPodesta et al.| 20061 

20071 iTessein et al. 20091 

Salem 

et al.j|:^UUt) Chen et aJ 

. 2011 bp, are very stable in time. 


T'hey appear at the maximum of the turbulent activity 
and persist throughout all the simulations, as the energy 
reservoir at lar ge scales is able to sustain and maintain 


the cascade. In [Letter ij we showed that such magnetic 
and velocity scaling are also combined with a spectral 
index of —2 for the residual en ergy, in agreement with 


observations in the solar wind dChen et al. 2Qjja). In¬ 


compressible MHD (|M ujjnr and Re 


duced MHD ( Boldyrev et al.||2(JIip only partially 


repro¬ 
duce such scaling. In our simulations, the 2D geometry 
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and the presence of compressibility may play a role in 
setting the different scaling. 

A clear transition in the spectra is observed at scales 
k_\_dp > 1, with a change in the spectral indices of all 
fields. In particular, the spectrum of the perpendicular 
magnetic fluctuations steepens at k± dp 2^ following a 
power law with spectral index ^ —3 for another decade. 
The location of the break does not show any significant 
dependence on the number of particles, the spatial res¬ 
olution and the resistivity adopted, provided that a suf¬ 
ficient number of grid points allows to cover approxima- 
tively a decade at sub-proton scales, i.e., that the scale 
at which resistive dissipation acts is sufficiently separated 
from the region of the break. The parallel component of 
the magnetic field, together with the density, follows a 
similar but slightly shallower slope with a spectral index 
of ^ —2.8, in very good agreement with obser vations 
(iChen et al. |2Q12 _ 2Q13b ) and o ther simulations (Howes 


et al. 2U11 "Fassot et al. 2014). As a result, magnetic 
huctnations tend to become isotropic at small scales, re¬ 
sulting in an increase of th e magnetic compressibility, as 


observed in the solar wind dSalem et al. [20121 [Podesta fc 
TenBarge|2012 Kiyani et al.|20l3| ). I’he spectrum of the 


perpendicular velocity huctuations quickly drops above 
k±dp ^ 1, without any clear power-law trend. The ob¬ 
servation of a spectral index of —2.8 has been ascribed 
to the effect of the electron Landau damping by previ¬ 
ous studies ( Howes et al.||2011 Passot et al.||2014 ), how¬ 
ever, this can not be the case in our simulations, where 
the electron kinetics is not taken into account. Alterna¬ 
tively, the presence of coherent structures, such as cur¬ 
rent shee ts, can produce a steepening of the energy spec¬ 
tra (e.g.,|Wan et al.||2012||Karimabadi et al.||2013). The 
increase of intermittency at small scales, observed in our 
simulations, seems to confirm this path towards the dis¬ 
sipation. We have to note, however, that a —8/3 power 
law for the magnetic energy and the density spectra (not 
far from the 2.8 found here) has been also interpreted as 
related to the dimensionality (ID or 2D) of the magnetic 
and the den sity intermittent structures, without invoking 
dissip ation (jBoldyrev & Perez||2012t |Meyrand & Galtier 

^M3). 

The spectrum of the electric fluctuations is highly dom¬ 
inated by its perpendicular component. It is strongly 
coupled to the spectrum of the perpendicular velocity 
fluctuations at fluid scales, then it decouples and flat¬ 
tens, exceeding the spectrum of the perpendicular mag¬ 
netic fluctuations and becoming dominant for k±dp ^ 2. 
At large scales, the only contribution comes from the 
MHD term in Eq. whose leading term is 

Ej_ (X u X B u± X Bq. (7) 

This corresponds to a power-law scaling 

Pe oc oc (8) 

which is observed in th e simulations, and is also consis¬ 
tent with observations ( Chen et al.||2Qlla ). In our case, 
the main contributions at sub-proton scales come from 
the Hall and the electron pressure gradient terms, since 
the spectrum of the velocity fluctuations is observed to 
drop exponentially at short wavelengths. The leading 
terms at these scales are, then. 


J X B — Vpe OC V(Ten + Bq • B\\ 


(9) 


Note that, in general, the sum of the two terms inside 
parentheses would not necessarily result in a power law 
for the electric field. However, since in our simulations 
both n and B\\ are observed to scale with the same power 
law - thanks to the strong coupling between the plasma 
and the magnetic compressibility - then the expected 
spectral index for the electric field is: 

Pe oc k‘^PB^^,n oc kj^'^ . (10) 

Although it is not possible to directly test this scaling 
for the electric field spectrum in Run A, individual terms 
in Eq. ^ follow well the prediction (see Eig. [^. More¬ 
over, we were able to show that when assuming Tg = 0 
(i.e., setting to zero the electron pressure gradient term 
in the Ohm’s law), then the electric field spectrum - 
hence dominated by the Hall term - follows a kj^'^ scal¬ 
ing in the sub-proton range (Eig. [^, being only very 
slightly affected by numerical effects at very small scales 
{k±dp > 10 ). 

As a result of the interaction of particles with the tur¬ 
bulent fluctuations and small scale structures, we observe 
an overall parallel and perpendicular heating with simi¬ 
lar rates, so that the temperature of the plasma remains 
globally nearly isotropic. This behavior can be achieved 
only if a high enough number of particles is employed, 
and the resistivity is properly set in order to assure an ac¬ 
curate conservation of the total energy and a clear power- 
law behavior for the spectrum of the magnetic fluctua¬ 
tions at all scales. The parallel temperature, T^n, is found 
to have a very robust evolution, being essentially inde¬ 
pendent of the resistivity, the number of particles, and 
the spatial resolution employed. On the contrary, the 
time evolution of Tp_\_ is strongly determined by both the 
resistivity and the number of ppc: if too few particles 
are employed, or if the resistivity is too low, the perpen¬ 
dicular heating can be largely overestimated/unphysical. 
Conversely, when a too strong value of the resistivity 
is implemented, the artificial damping of fluctuations at 
ion scales can produce a strong reduction of the perpen¬ 
dicular heating, thus generating an equally unphysical 
preferential parallel heating. This proves that no firm 
conclusions can be drawn about the perpendicular heat¬ 
ing by turbulence in hybrid simulations, unless a careful 
and empirically fine-tuned choice of all parameters has 
been taken. 

Note, however, that the fact that we do not observe a 
global preferential heating does not imply the absence 
of signatures of localized preferential deformations of 
the particle distribution functions, as suggested by the 
bottom right panel of Eig. where strong temperature 
anisotropies ranging from 0.5 to 1.8 are observed. They 
seem to be concentrated in regions with stronger coherent 
structures, identified by the presence of current sheets 
and a significant level of vorticity. These results are 


in agreement with previous 

works based on the Vlasov- 

hybrid approximation (e.g., 

Servidio et al.||2012l 

Perrone 

et al. 2013 Valentini et al. 

2014| bervidio et aJ 

f: Mir 


As the overall heating is rather weak, slow, and nearly 
isotropic, we can infer that the local formation of large 
proton temperature anisotropies is likely due to energy 
exchanges between the parallel and perpendicular direc¬ 
tions, and/or to the spatial transport, rather than due 
to the heating. 
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Solar wind observations show a certain variability of 
the spectral properties. In particular, the position of the 
break at ion scales and the shape of the magnetic field 
spectrum around i t seems to depend o n the power of mag¬ 
netic fluctuations (Bruno et al.||2Q14|) and on the plasma 
beta ([Alexandrova et ah ' 2 UU8a[ ~ |Chandran et al. 2009 


Chen et al.|2014p . investigating such a dependence, by 


exploring the parameter space of the level of fluctuations 
and the plasma beta, will be the subject of a fortcoming 
paper. 

Three-dimensional simulations would be fundamental 
to further improve the present study and overcome its 
limitations, allowing for a more realistic description of 
the turbulent cascade. 
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